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Abstract. We investigate the thermal, magnetic and ro- 
tational evolution of isolated neutron stars assuming that 
the dipolar magnetic field is confined to the crust. Our 
treatment, for the first time, uses a fully general relativis- 
tic formalism not only for the thermal but also for the 
magnetic part, and includes partial general relativistic ef- 
fects in the rotational part. Due to the fact that the com- 
bined evolution depends crucially upon the compactness 
of the star, three different equations of state have been 
employed in the calculations. In the absence of general 
relativistic effects, while upon increasing compactness a 
decrease of the crust thickness takes place leading into an 
accelerating field decay, the inclusion of general relativistic 
effects intend to "decelerate this acceleration" . As a con- 
sequence we find that within the crustal field hypothesis, 
a given equation of state is compatible with the observed 
periods P and period derivative P provided the initial field 
strength and current location as well as the magnitude of 
the impurity content are constrained appropriately. 

Finally, we access the fiexibility of the soft, medium 
and stiff classes of equations of state as candidates in de- 
scribing the state of the matter in the neutron star inte- 
riors. The comparison of our model calculations with ob- 
servations, together with the consideration of independent 
information about neutron star evolution, suggests that a 
not too soft equation of state describes neutron star interi- 
ors and its cooling proceeds along the 'standard' scenario. 
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1. Introduction 

The question whether and, if so, how the magnetic field 
of isolated neutron stars (NSs) decays is a controver- 
sial issue and a subject of hot scientific debates. The 
observed rotational periods P and time derivative P of 
about 700 pulsars (PSRs) (see PSR catalogue, Taylor et 
al. 1993) and studies of their inferred surface magnetic 



field strength versus their active age (ra) provide evidence 
that the magnetic fields of NSs is subject to decay. This 
evidence is rather strong for the old population of PSRs, 
i.e. Ta ^ 100 Myrs, while for the younger population evi- 
dence for magnetic field decay is much weaker. As investi- 
gated by population synthesis methods (Bhattacharya et 
al. 1992, Hartman et al. 1997), observations and mod- 
els are in harmony provided one accepts the hypothe- 
sis that NS magnetic field decays very little during the 
first Ta ^ 10 Myrs. The NS magnetic field decay depends 
strongly on where is the field located in the NS, what is its 
structure and strength. It is also related to the equation of 
state (EOS) of NS matter and the conductive properties of 
dense matter which are moreover affected by its thermal 
history too. In case the field is of fossil origin, or alter- 
natively has been generated via a dynamo action during 
the proto-NS phase (Thompson & Duncan 1993), it is ex- 
pected to thread most of the star. Early estimates of the 
electric conductivity of NS matter lead to the conclusion 
that magnetic field of NSs will not be dissipated during 
a Hubble time (Baym, Pethick & Pines 1969). However, 
more recent investigations regarding the conductive prop- 
erties of nuclear matter in the presence of a strong B-field 
yield conductivities leading to much shorter field decay 
times (see e.g. Haensel, Urpin & Yakovlev 1990, Goldre- 
ich & Reissenegger 1992, Urpin & Shalybkov 1995, Shaly- 
bkov & Urpin 1997). Besides the conductive properties of 
nuclear matter, another distinct mechanism leading to a 
B-decay process is based on the idea of magnetic flux ex- 
pulsion from NS core driven by rotation or/and buoyancy 
(Muslimov & Tsygan 1985; Srinivasan et al 1990). Once 
the expelled field reaches the star's crust it subsequently 
suffers Ohmic decay. Ding et al. (1993) considered in de- 
tailed this process and have found typical decay times of 
the order of 100 Myrs. Recent investigations, which take 
into account the effect of the NS crust onto the process 
of flux expulsion (Konenkov & Geppert 2000) estimated 
even longer decay times for a field anchored in fluxoids in 
the superfluid NS core. 
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Although the issue of whether the field is penetrating 
the entire star or part of it is an open one, there exist, 
however, good reasons to believe that the NS magnetic 
field is maintained by currents in the crust. One possi- 
bility is the generation of the field via thermomagnetic 
effects during the first years of the NS life (Blandford ct 
al. 1983, Urpin, Levshakov & Yakovlev 1986, Wiebicke & 
Geppert 1996). The most recent and detailed investiga- 
tion of the magnetic and spin evolution of isolated NSs 
with a crustal dipolar field has been performed by Urpin 
& Konenkov (1997, UK97 thereafter). They found a good 
agreement with observational data provided the EOS is 
not too soft, the initial surface magnetic field strength lies 
in the range of 10^^ to 3 x 10^'^ G and is initially con- 
fined to densities of 10^^ to 10^"^ g cm^"^ within the crust. 
Miralles, Urpin & Konenkov (1998) considered the effect 
of a crustal field decay onto the thermal evolution of a 
NS and they have concluded that a considerable amount 
of heating takes place in the the crust after about 3 to 
10 Myrs , a period during which the NS has cooled down 
to ^ 10^ K. Consequently, Joule heating can maintain a 
warm (« 5 ■ 10^ K) NS surface for a period of hundreds of 
Myrs. 

All of the above described works ignore the curvature 
of spacetime. To date there exist no calculation treating 
self consistently the thermal and magnetic evolution of 
NSs which incorporate General Relativistic (GR) effects. 
As it has been argued elsewhere by the present authors 
(Geppert, Page & Zannias 2000, referred to as GPZOO 
hereafter), relativistic effects on the B-field evolution must 
be included in detailed investigations. A previous attempt 
to incorporate GR effects has been presented by Sengupta 
(1997, 1998). However this author failed to take into ac- 
count the proper boundary condition associated with rel- 
ativistic treatment and furthermore his formalism applies 
only to a Schwarzschild geometry. Moreover, he claims 
that the decay rate of the field is decreased by several or- 
ders of magnitude, a conclusion being in variance with the 
one obtain by GPZOO and the present work. 

In the present paper we investigate in details and in 
a self consistent manner, GR effects at first on both: the 
thermal and magnetic evolution of NSs. We consider three 
different EOSs and via numerical integration of Einstein's 
equations, neutron star models are constructed character- 
ized by the compactness ratio varying in a large range. 
We solve simultaneously the relevant evolution equations 
and thus our approach naturally reveals the mutual de- 
pendencies of EOS, mass-to-radius relation, initial field 
structure and strength, and thermal history of NSs. In 
particularly, as far as the cooling is concerned we consider 
two scenarios: the 'standard' slow neutrino emission sce- 
nario and also the so called enhanced neutrino emission 
which results in a much lower temperatures in young NSs. 
In order to avoid dealing with the uncertainties related to 
the behavior of the field within the superconducting core, 
as a first step in the present paper, we consider only mag- 



netic fields not penetrating the NS core. Thus, the present 
analysis is within the framework of the crustal magnetic 
field hypothesis. We use our results from the GR treatment 
of thermal and magnetic field evolution implemented by 
a semi-relativistic treatment of the spin evolution to con- 
front the crustal field hypothesis to observational data. 

The paper is organized as follows: In the subsequent 
section we remind the reader of the GR formulation of the 
equations of stellar structure, the heat transport and con- 
servation equations as well as the induction equation on 
a static spherically symmetric background geometry. We 
also present the evolution equations for an axisymmetric 
dipolar magnetic field as well as relativistic expressions for 
the Joule heating. Section 3 describes the microphysics 
used in our models. In section 4 we present the results 
of our model calculations and section 5 is devoted to the 
discussion of the results. 

2. General Relativistic formalism of static 
spherical stars 

The spacetime geometry will be assumed to be spheri- 
cally symmetric, which means that our results may not 
be accurate for fast rotating neutron stars. Employing the 
familiar Schwarzschild coordinates (i,r, 0, (/>), the interior 
and exterior spacetime geometry takes the form (Misner 
et al. 1973; Wald 1984) 



1 - IGmjc^r 



+ r^dn^. 



(1) 



The radial proper length is thus dl = drj^X — 2Gml(?r 
and the proper time dr = dt. 

Einstein's equations coupled to a perfect fiuid energy- 
momentum tensor give us the standard equations (Misner 
et al. 1973; Wald 1984): 



dm 2 
— — = 47rr p, 
dr 



(2) 



which determines the so called mass function m ~ m(r) 



dr c4r2(l -2G'm/c2r)' 



(3) 



for the 'gravitational potential' <& — $(r) and the TOV 
equation of hydrostatic equilibrium 

dr dr 



{p + P/^){Gm + A-KGr'^P/c^) 
r2(l - IGmjc^r) ' 



(4) 



Regularity of the geometry at r = implies the inner 
boundary condition for Eq. ^ 



m(r = 0) = 

while for Eq. |^ the central pressure 
P(r = 0) = 



(5) 
(6) 
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is specified, through the EOS of the form P = P{p), with 
the central density pc treated as a free parameter. Due 
to the linear nature of Eq.(3), $ can be scaled so that it 
always can be arranged to fulfill: 



1 



2GM 



(7) 



Once the interior spacetime geometry has been so specified 
it is joined smoothly across the "surface" of the star to an 
exterior Schwarzschild field characterized by M = m{R). 

It should be mentioned that the stellar surface in our 
computation is fixed by 



which relates the temperature at the outer boundary (de- 
fined more precisely further bellow,), Tf,, to the luminosity, 
it,, in this layer. The location of this outer boundary layer 
is chosen such that Lt be equal to the total photon lumi- 
nosity of the star, = L{r = R), which in turn is related 
to the effective temperature Te by i* = AnR^asBT^- We 
can thus write Eq. |l| as Tf, = Tb(Te) and this 'Tf, - 
relationship' is discussed in Sect. 



2.4 



The opacity is related to the total (thermal) conduc- 
tivity by: 



A = 



■ Hp = Pb) (8) 

10^" g cm"-^. This guarantees that the EOS 
is temperature independent. The layers at densities below 



R = Rs 



where pi, 



Znp 



(14) 



If one neglects the red-shift e* and defines the energy fiux 
F as L/{A'Kr'^) one can write Eq. p2| as 



pb, called the envelope, are treated separately (see § 2.4). p = —\ . SJT 



(15) 



2.1. Thermal evolution equations 

Besides the above equations of stellar structure we shall 
need the equations describing the thermal evolution of the 
star. At the temperatures we are interested in, the neutri- 
nos have a mean free path much larger than the radius of 
the star (Shapiro & Teukolsky 1983) and thus leave the 
star once they are produced. Energy balance arguments 
(see for instance Thorne (1966) then imply 



d(LO 
dr 



Anr'^ne^ 



- IGvajc? 



(9) 



where L is the internal luminosity, e the internal energy 
per baryon, q^^qh is the neutrino emissivity and heating 
rate per baryon while n stands for the baryon number 
density. The corresponding inner boundary condition for 
L is 

L(r==0) = (10) 
The time derivative of e can be written in the form 



de de dT 
dt ~ dr ' ~dt 



dT 



(11) 



through the specific heat at constant volume (which, 
for degenerate matter, is the same as the specific heat at 
constant pressure cp). 

The energy transport equation is : 



d{Te^ 
dr 



up 



Le'' 



(12) 



where T is the local temperature, asB the Stefan- 
Boltzmann constant and k, the 'opacity'. (Notice that 
within the relativistic framework an 'isothermal' config- 
uration is defined by e"^ • T = constant instead of T = 
constant.) The associated boundary condition is 



Tb — Tb{Lb 



(13) 



where V is the radial gradient calculated with the proper 
length, which is the usual form of the heat conduction 
equation. The total conductivity is the sum of the electron 
and photon conductivities 



A — Ae -|- A^ 



(16) 



since these two processes of heat conduction work inde- 
pendently and in parallel. 

We will present our results of thermal evolution by 
using the 'effective temperature at infinity' T^ = Tg-e**^^^ 
related to the 'luminosity at infinity' e 
through the 'radius at infinity' R°° = R ■ e^*'^^^ by 



(17) 



These three quantities 'at infinity' are, in principle, mea- 
surable and, in particular, R°° would be the areal radius of 
the star that an observer 'at infinity' would measure with 
an extremely high angular resolution instrument (Page 
1995). 

The solution of the complete set of equations of stellar 
structure, Eqs. (2-7) and thermal evolution, Eqs. (9-13), 
requires knowledge of the equation of state, the opacity 
K, the neutrino emissivity and also the specific heat 
Cy. We shall devote Sect. ^ to the detailed specification of 
those variables. 



2.2. Magnetic evolution equations 

Besides the equations of stellar structure and thermal evo- 
lution we shall also need the equations describing the mag- 
netic field evolution. In this paper consideration will be 
restricted to dipolar (poloidal) magnetic fields. The GR 
formulation of the evolution equation of such fields has 
been discussed in detail by GPZOO while for a more gen- 
eral set up see Radler et al. (2000). As has been shown in 
the first reference, such a field can be expressed in terms of 



4 



Page, Geppert, & Zannias: Isolated Neutron Stars with Crustal Fields 



a relativistic generalization of the familiar Stoke's stream 
function by: 



B'{t,r,9) = ^^^^ COS0, 



r \ err J 



1 2GmY dF{t,r) 

dr 



sm( 



(18) 
(19) 



while the relevant induction equation (see Appendix) 
yields the following equation for the relativistic Stoke's 
function 



Ana 



dF 
'dt 



= 1- 



2Gm\ d^F 
c?r I dr-^ 



+ 



1 dF 

J.2 Qj. 



2Gm A-kG 



tF. 



(20) 



The appropriate boundary conditions, as r — > 0, is the 
same as in the flat space case: a regular field at the star's 
center requires ^'^^■f'^ to be finite. The outer boundary 
condition, however, differs from that valid in the flat space, 
and its GR form is as follows (see: GPZOO): 



R 



dF{t,r) 



dr 



Giy)Fit,R) 



where: 



Giy) 



2yln(l-r 



2jy-l 



^ y"^ ln(l — y ^) 



y 



(21) 



with y = R/Rs (22) 



(i?S = 2GM/c^ being the star's Schwarzschild radius). 
Since G(y) < (in particular, in the flat space-time case, 
G{oo) — —1) the boundary condition forces a bending of 
F in the upper layers. 

As an initial profile for the Stoke function we use the 
same formula as UK97, for later comparison with the work 
of these authors: 

i^(r,0) = Boi?2 (l-rVr^)/(l-i?Vr§) at r > ro (23) 



F(r,0) = 



at r < Tq 



Notice that this initial F{r,0) does not satisfy the outer 
boundary condition, but will immediately be forced to do 
it at the first numerical time step. There will thus be a 
rapid relaxation of F in its early evolution due to the en- 
forcement of the boundary condition and the propagation 
of the resulting curvature of F toward higher densities. 

2.3. Joule heating 

Due to the finite conductivity cr, magnetic energy is dis- 
sipated into heat (= Joule heating). The heat production 
per unit (proper) time and unit (proper) volume is given 

by 



Qh = n ■ qii = — 



(24) 



where j is the current. As it is shown in the Appendix, for 
a dipole poloidal magnetic field B, j is given by j = i^e^ 
with 



c sm( 



J0 



47r 



1 



2Gm'Y d 
c^r J dr 



1 



2Gm 



WF 
dr 



2F\ 



(25) 



It is seen immediately that in the limit e* = 1 and m = 
the above expression reduces to its flat space-time coun- 
terpart (Miralles et al. 1998). Since our numerical calcu- 
lations assume spherical symmetry we use, in Eq. |2^, a 
spherical average of sin^ 9 (< sin^ 9 > — |) obtaining 



< Qh >-- 



1 1 



247r2 , 



1 - 



2GmV- 
(?r j 



d_ 

dr 



2Gr. 



WF 
dr 



2.4- Outer boundary: magnetized envelopes 

The layers at densities below pf, = 10^*^ g cm~'^ are defined 
as the envelope, and extend up to the atmosphere, where 
the photosphere is located, while by interior we mean the 
whole star where p > pb- The presence of the magnetic 
field affects strongly the heat transport in the envelope 
(but not in the deeper layer where p > 10^° g cm""^, 
which motivates the above choice of pb), and results in a 
non uniform distribution of the surface temperature. The 
corresponding Tf, — relationships' have been calculated 
by Page & Sarmiento (1996) for dipolar (and quadrupo- 
lar) fields. We use these results which hence give us a field 
dependent Tf, — Tg relationship' that adjust itself to the 
evolution of the magnetic field. Notice however that, for 
fields much stronger than 10^^ G, this relationship is not 
reliable when is much lower than 10^ K, and it is most 
probably very inaccurate when Tf, is below 10^ K. 

Our outer boundary condition assumes that the en- 
velope is made of catalyzed matter. If an upper layer of 
light elements were present, the heat transport is strongly 
enhanced when no magnetic field is present (Potekhin, 
Ghabrier, & Yakovlev 1997), but there is, to date, no pub- 
lished model of magnetized envelopes with light elements. 

One should finally emphasize that when the thermal 
evolution is controlled by the Joule heating the star's lu- 
minosity, and Te, is given by the heating rate and is inde- 
pendent of the outer boundary condition as discussed at 
the end of '• 



4.2 



This closes the system of equations and boundary con- 
ditions to be solved. 

2.5. Numerical method 

The thermal evolution equations are solved by a Henyey- 
type code (e.g.. Page 1989) while the induction equation 
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for the Stoke function is solved with a Crank-Nicholson 
method (Press et al., 1986). The whole set of equations 
for the thermo-magnetic evolution should be solved si- 
multaneously at each time step but we have decided to 
solve the thermal and then the magnetic equations alter- 
natively, i.e., the thermal equations are solved at a given 
time step using the field of the previous step (which ap- 
pears in the Joule heating term) and once the new tem- 
perature is obtained the induction equation is solved to 
obtain the new magnetic field. This method is much faster 
than a full simultaneous solution and gives results which, 
as we have verified explicitly in a few case, are practically 
indistinguishable from the full simultaneous solution. 

3. Input microphysics 

3.1. The equation of state 

The first ingredient needed to build NS models is the equa- 
tion of state (EOS). In principle the EOS should give us 
not only the relationship between pressure and density, 
i.e., P = P(yo), but also the chemical composition of mat- 
ter. 

We separate the crust from the core at the density 
p = Per = 1.6 X 10^'' g cm""^ (Lorenz, Ravenhall & Pethick 
1993) and use the EOSs of Negele & Vautherin (1973) for 
the inner crust, at p > pdrip = 4.4 x 10^^ g cm^'^, and 
Haensel, Zdunik & Dobaczewski (1989) for the outer crust, 
i.e., we assume that the chemical composition is that of 
cold catalyzed matter. The EOS, and its associated chem- 
ical composition, in the crust is well determined under 
the assumption that matter is in its (catalyzed) ground 
state. There is however still the possibility that a strong 
phase of hypercritical accretion occurred after the super- 
nova explosion, which may, or may not, alter the chemical 
composition of the crust. We will not consider this possi- 
bility here but only mention that it does have an enormous 
effect on the magnetic field and its subsequent evolution 
(Geppert, Page & Zannias 1999) 

In the core, the EOS is relatively well constrained up 
to densities around 2 - 3 x pcr while its behavior at higher 
densities is still a mystery (Prakash 1998). We will thus 
consider three different cases which hopefully illustrate the 
whole range of possibilities. We take as a 'Medium EOS' 
the one calculated by Wiringa, Fiks & Fabrocini (1988), 
using their model called avl4+UVII, a 'Stiff EOS' from 
Pandharipande, Pines & Smith (1976) and a 'Soft EOS' 
from Pandharipande (1971). Notice that the Stiff EOS is 
based on the presence of a lattice of neutrons in the in- 
ner core and is not anymore considered as realistic but we 
still use it since it has been used by many authors and 
represents the case of extreme stiffness. Moreover, none of 
the consequences of the presence of such a lattice phase is 
taken into account in our calculations, e.g., in C„, and 
the transport coefficients: we will assume that neutrons 
and protons form a quantum liquid and also boldly con- 
sider them as superfluid and superconductor. Our opinion 
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EOS 


Radius 

[km] 


Crust 
thickness 

[km] 


Crust 
mass 
[Me] 


Central 
density 

[10^5 
g/cm^] 


Moment 
of inertia 
[10^=^ 
g cm^] 


Stiff 


14.94 


1.52 


0.069 


0.49 


1.83 


Medium 


10.48 


0.72 


0.022 


1.17 


1.08 


Soft 


7.1 


0.22 


0.003 


5.59 


0.72 



Table 1. Properties of the 1.4 Mq neutron star models 



is that this EOS should be abandoned but we consider it 
for comparison with previous works of other authors. The 
Soft EOS on the other hand, despite of its age, is still rep- 
resentative of modern soft EOSs and includes hyperons. 

We will consider models of 1.4 Mq stars whose overall 
properties are listed in Table H. 



3.2. Neutrino processes 

Neutrino emission drives the cooling as long as the internal 
temperature is higher than about 10® K. For processes in 
the crust we consider the two dominant ones which are 
the plasmon process and the electron-ion bremsstrahlung 
(Page 1989). The former is only relevant during the first 
few years of the life of the NS but is very strong and brings 
down the crust temperature to about 10^ K, relaxing it 
from the arbitrary initial conditions. The latter has only 
a very small effect, mostly when the surface temperature 
is around 10^ K, i.e., for young stars. 

The crucial neutrino emission processes occur in the 
core and we consider two scenarios (see Page 1998 for more 
details). In the 'Standard Cooling', or 'slow cooling', sce- 
nario we include the modified Urea processes and their 
associated, and weaker, bremsstrahlung processes follow- 
ing Yakovlev & Levenfish (1995). This scenario applies 
when the NS core contains only neutrons, protons, as well 
as electrons and muons which maintain charge neutral- 
ity, and the proton fraction is low enough that the di- 
rect Urea process is forbidden by momentum conserva- 
tion. In the 'Enhanced Cooling', or 'fast cooling', cases 
we add a strong neutrino emission at densities larger to 
Pfast = 4 X 10^^ g cm^'^ with a rate 

= lO^s ip/perr/' n erg g-1 s-1 (27) 

where Tg = T/10^ K. This rate is representative of many 
of the possible enhanced ones as the direct Urea from nu- 
cleons or hyperons, but is stronger than what produced 
by a pion or kaon condensate (Prakash 1998). For, com- 
parison the inefficient modified Urea process gives approx- 
imately e^URCA _ ^q21 rj^l g-1 

In our three model 1.4 Mq stars, the masses of the 
inner cores where the fast neutrino emission is allowed 
for the 'Enhanced Cooling' cases are, 0.11, 1.29 and 1.39 
Mq for the Soft, Medium and Stiff EOS, respectively. 
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10 



Density [g cm"^ ] 

Fig. 1. Electrical conductivities in the crust. The dashed 
lines show crimp, which is temperature independent, for 
three different values of the impurity content Qimp- The 
continuous lines show a for electron-phonon scattering, in 
the solid phase, and electron-ion scattering, in the liquid 
phase, at four different temperatures. Notice that these 
four curves all correspond to the solid phase except for the 
10^ K one where the low density part is still in the liquid 
phase: the dotted extension show the value a would have 
if the matter were in the solid phase and the short tripped 
segment indicates the region where a is interpolated be- 
tween the electron-phonon and the electron- ion values. Fi- 
nally, the increase of a in this 10^ K case at densities just 
above the transition density, i.e., temperatures just be- 
low the melting temperature, is due to the Debye- Waller 
factor (Itoh et al 1984). For comparison with the cooling 
curves, notice that interior temperatures of 10^, 10^, 10^, 
and 10^ K roughly correspond to effective temperatures 
around 10^, 3 x 10^, 10^, and 3 x 10^, respectively. 
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Fig. 2. Pairing critical temperatures adopted in this work, 
for neutrons in the ^Sq state (from Ainsworth, Wambach 
& Pine 1989) and state (from Takatsuka 1972) and 
protons in the ^Sq state (from Baldo et al. 1992). The 
dashed lines show the central densities of our three model 
1.4 Mq stars with the soft, medium and stiff EOS. 



T = To_i for electron-ion scattering by Itoh et al (1983). 
In the solid phase t is given by 

1 1 1 



-ph 



(30) 



where we use the results of Itoh et al (1984) for the 
electron-phonon scattering rc_ph and of Yakovlev & Urpin 
(1980) for the electron-impurity scattering Tc-imp. In gen- 
eral, Te-i is a function of p, A and Z, and Tc_ph depends 
also on the temperature T while rc_i,np depends on p, Z 
and the impurity concentration Qimp- We show in Fig. |^ 
the value of a for typical values of T and Qimp- 



3.3. Electrical conductivities 

In general, the electrical conductivity tr can be expressed 
in terms of the electron relaxation time r as 



a = 



where 



^^e/c^ = me [1 + 1.018 (pe ZjAfl^X^I"^ 



(28) 



(29) 



is the electron effective mass and Ue the electron number 
density (/Xe being the electron chemical potential, p^ = 
p/10^ g cm~'^, Z and A the charge and mass number of 
the ions). In the liquid phase we use the calculation of 



3.4. Pairing 

The occurrence of pairing, neutron superfluidity and pro- 
ton superconductivity, strongly affects both the thermal 
and magnetic evolution of neutron stars. The thermal ef- 
fects are very strong during the neutrino cooling phase, 
which last about 10^ to almost 10^ yrs depending on the 
model, and the subsequent photon cooling phase (see, e.g.. 
Page 1998 for a review). As a result, pairing will in a large 
part determine the time at which Joule heating starts to 
control the thermal evolution and then during this Joule 
heating phase the effect of pairing becomes negligible. We 
treat the suppressive effect of pairing on Cy and ei, ac- 
cording to the treatment of Levenfish & Yakovlev (1994a, 
1994b) and Yakovlev & Levenfish (1995). 

The superconductive phase in the core has the effect of 
producing an impenetrable barrier for the magnetic field. 
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Model 


Bo [G] 


Po [gm/cm^] 


Qimp 


1 


10" 


10^2 


0.1 


2 


10^2 


10" 


0.01 


3 


10" 


10^^ 


0.001 



Table 2. Initial magnetic field configurations and impu- 
rity concentrations in our three classes of models. 



which is initially confined to the crust in the models of 
the present work. This guarantees the confinement of the 
magnetic field, and the currents, to the crust since the su- 
perconductive phase transition happens, at the crust-core 
boundary, well before the magnetic field had time to dif- 
fuse to this layer. For definiteness we plot in Fig. ^ the 
pairing critical temperatures that we adopt in this paper. 
Notice that the values we adopt for ^Sq pairing of both 
neutrons and protons are typical of modern calculations. 
In the case of neutron pairing we explicitly adopt val- 
ues which vanish a densities above 10^^ g cm~^ to ensure 
that pairing will not suppress the strong neutrino emis- 
sion in our Enhanced Cooling scenarios for the Soft and 
Medium EOSs. 

In the case of the Stiff EOS any published calculation 
of neutron or proton pairing shows a non vanishing value 
of Tc at the density in the center of the star, given the low 
value of this central density. To avoid suppression of the 
neutrino emission by pairing we will assume that the fast 
neutrino emission is not affected by neutron and proton 
pairing in this case of Stiff EOS in the Enhanced Cooling 
scenario . 



4. Results 

We will present here our results for the thermal, mag- 
netic and rotational evolution of isolated neutron stars 
considering 1.4 Mq stars built with the three EOSs de- 
scribed above within both the 'standard' and the 'fast' 
cooling scenarios. In order to investigate the influence of 
the (a priori unknown) initial structure and strength of 
the magnetic field onto NS's evolution, we considered for 
each EOS and cooling scenario three classes of qualita- 
tively different field models, existent at the beginning of 
the evolution. They are characterized by the initial surface 
field strength Bq , depth of penetration of the current po 
and impurity concentration Q as listed in Table 0. 



A detailed example 

We first discuss here an example of the internal evolu- 
tion of the Stoke function, magnetic field, currents and 
local heating rate, as shown in Fig. |^, which will help for 
the general discussion presented below. We choose the 1.4 




p [gm/cnnS] 



Fig. 3. Stoke function, magnetic field, B, currents, j and 



Joule heating rate, 
central panel of Fig. 



in the crust, for model 2 of the 
See text, §4.1, for details. 
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Age [yrs] 



Fig. 4. Evolution of the internal temperature (lower 
panel) and surface magnetic field for the model described 
in Fig. 1^, continuous line. The two models in dashed and 
dotted lines show the corresponding evolution for impurity 
contents Qimp = 0.001 and 0.1, respectively (the continu- 
ous line has Qimp — 0.01). Notice that in these cases the 
heating is so small that the three models follow the same 
thermal evolution track. 



Mq star with the medium EOS and 'standard' cooling; the 
initial surface field strength is 10^^ G and the currents are 
initially located at po = 10^'^ g cm~^. The impurity con- 
tent for (Tinip is Qimp = 0.01. This is the model 2 of the 
central panel of Fig. ]| but it is reproduced in Fig. ^ along 
with two similar models with different impurity contents 
(upper panel) and the evolution of the internal temper- 
ature (lower panel). GR effects are included but will be 
discussed in the next subsections. 

The four chosen times, labeled as a, b, c & d are 
marked in the upper panel of Fig. ^ and correspond, re- 
spectively, to a) the initial field decay during the period 
when the neutron star is still hot and a is temperature 
dependent, being dominated by electron-phonon scatter- 
ing, b) the plateau where the cooling lead to an enor- 
mous increase of a, and consequently a stagnation of the 
field decay, which will eventually be controlled by impu- 
rity scattering, c) second phase of decay of the field when 
time becomes comparable to the impurity scattering decay 
time scale and, finally, d) the late exponential decay. 

The four panels of Fig. || directly illustrate the dif- 
fusion of the field toward regions of higher conduc- 
tivity as time runs. Once a becomes T independent 
the diffusion equation formulates an eigenvalue problem 
whose solution can be formally written as F{p,t) = 
'3xp(-t/r„)a„F„(/9), in terms of eigenmodes F„ 
with decay times t„, and expansion coefficients a„. Since 
the n-th mode, F„{p), has n nodes, in the crust, initially 
a very large number of modes must contributes signifi- 



cantly, i.e., with large positive and negative coefficients 
a„, to produce a mutual cancellation resulting in a vanish- 
ing F{p, t) in the high density region. The diffusion of the 
stoke function into the high density region is simply due to 
the faster decay of the modes with n > 1 compared to the 
nodeless fundamental mode. When the fundamental mode 
is dominating, i.e., when the stoke function is non zero in 
the whole crust (it has reached the crust-core boundary 
through diffusion), the field evolution is a power-law like 
decay, phase c (Urpin, Chanmugan & Sang 1994). Finally, 
when only the fundamental mode n = 1 is left the decay 
becomes purely exponential, phase d. 

Notice that the field strength inside the crust is at least 
one order of magnitude higher than at the surface: this is 
due to a very large Bg forced by the presence of the deriva- 
tive dF/dr in Eq. |l^, while at the surface the boundary 
condition ensures that Br and Bg are comparable. 

With respect to the currents, noticeable are the nega- 
tive currents in the layers where the field is growing which 
are due to induction, i.e., Lenz law. Once the field has 
reached the crust-core boundary the currents are positive 
in the whole crust and the only negative currents left are 
the supercurrents induced in the skin layer of the proton 
superconductor. Instead of imposing an ad- hoc boundary 
condition at the crust-core interface to simulate the effect 
of the proton superconductor we have preferred to keep 
the central boundary condition and introduce an enor- 
mous value for a, 10^'^'^ s^^, once protons become super- 
conductor, i.e., when T < Tc- This allows us to see explic- 
itly the induced supercurrents which, however, because of 
the finite radial resolution of the numerical scheme, are 
located in the whole last zone of the crust instead of the 
physical skin layer of the superconductor which is a few 
tens of fermis thick. We have checked explicitly that an 
ad-hoc crust-core boundary condition for F gives the same 
results as our boundary condition with enormous a. 

As a last remark, from Fig. |^, and comparing with 
Fig. |l|, we see that impurity scattering starts to dominate 
a after the field has reached the plateau (phase b): this 
stagnation value of the field in independent of Qimp and 
is only a result of the enormous increase of a due to the 
cooling. After this, the length of the plateau (phase b) 
is controlled by Qimp, higher impurity contents leading 
naturally to an earlier onset of the second phase of field 
decay (phase c). 

4-2. Magneto-thermal evolution within the 'standard' 
cooling scenario 

In Figs. H and ^ we present the modeled evolution of the, 
more or less, observable quantities, i.e. the surface tem- 
perature (measured by a far distant observer) and the 
surface magnetic field. While the surface temperatures, 
or, at least, its upper limits, of isolated NSs can be in- 
ferred from X-ray spectra, the magnetic field of isolated 
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Fig. 5. Thermal and magnetic evolution within the 'standard' cooling scenario. Evolutionary curves are labeled fol- 
lowing notations of table ||. See text for details. 



pulsars is mostly estimated by the precise measurement of 
the rotational period and its time derivative. 

Within the 'standard' cooling scenario the surface tem- 
perature stays quite high (~ 10^ K) during the neutrino 
dominated cooling era, which lasts about ^ 10^ yrs (for 
our stiff EOS models) to almost ~ lO'^ yrs (soft EOS mod- 
els). Later on the cooling is driven by photon emission 
from NS's surface which appears in the cooling curves by 
a strong increase in the slope. During the neutrino cool- 
ing era, the surface temperature drops down by a factor of 
about 2 (stiff EOS) to 4 (soft EOS), i.e. the approximately 
isothermal crust has at the end of the neutrino cooling 
epoch a temperature of about 5 • 10'' (soft EOS) to 10*K 
(stiff EOS). As seen from Fig. for such and higher tem- 
peratures, the electrical conductivity in the crust is deter- 
mined by electron-phonon collisions except at the highest 
densities. With the following temperature drop during the 
photon cooling era, the electron-phonon relaxation time 
Te_ph increase dramatically and thus impurity scattering 
will begin to control a (see Eq. 30 ) , at a temperature, and 
thus an age, depending on the impurity concentration. 
From that stage on, a becomes temperature independent. 



Notice that the differences in the thermal evolution 
for the different EOSs, before Joule heating becomes effi- 
cient, are mostly due to the differences in the fraction of 
the core which is paired. A larger paired region implies a 
lower neutrino emission, and thus a higher temperature 
during the neutrino cooling era, and also a lower specific 
heat, and thus an earlier transition to the photon cooling 
era and a faster temperature drop during that era (see, 
e.g.. Page 1998). We have taken here the choice of well 
defined density dependences of Tc, for both the neutrons 
and protons, independently of the EOS, as shown in Fig. |[ 
Different density dependences of the T^s would obviously 
give different results. For example, assuming high values 
of Tc down to center of the star for the soft EOS would 
results in a cooling history practically indistinguishable 
form the cooling history of the stiff EOS model. Given 
the present uncertainty on the value and density depen- 
dence of Tc for neutron pairing (Baldo et al. 1998) 
any choice has, unfortunately, some arbitrariness. Any ef- 
fect of the resulting cooling history on the field evolution 
should thus be specifically formulated in terms of cooling 
history and not in terms of the stiffness of the EOS. 
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Fig. 6. Cooling and field evolution within the 'fast' cooling scenario. Evolutionary curves are labeled following notations 
of table 0. See text for details. 



We have, very roughly, for the field decay time-scale, 
Tdecay oc P , I being a typical length scale of the crustal field 
structure. This immediately implies that Tdecay increases 
when pq increases and also when the stiffness of the EOS 
is increased, since both increase I. Moreover, since also 
Tdecay oc CT, a higher po locates the currents in a region 
of higher a and increases more Tdecay Thus, the models 1 
give fast decay, models 2 intermediate decay and models 
3 slow decay, as is clear form Fig. ||. 

The cooling influences the field decay by rising a till 
it becomes temperature independent when dominated by 
impurity scattering. This happens during the photon cool- 
ing era as mentioned above, and happens earlier for stiffer 
EOSs given our choice of EOS independent Tc 's. Conse- 
quently, we obtain that a stiffer EOS results in a slower 
field decay because of its cooling behavior and also be- 
cause of the larger length scale I. Our choice of Tc thus 
maximalize the effect of the EOS's stiffness on the field 
decay. 

We now turn our discussion to the analysis of GR ef- 
fects upon magnetic field evolution. For each model we 
simulated the field evolution with and without GR effects. 



but note that GR effects on the star structure and cooling 
have always included. (Had we for instance, have turned 
off the GR effects on the star structure, the resulting mod- 
els would make no sense at all since the differences in the 
size, the central density etc, of the model would tender 
them unrealistic.) On the other hand, GR effects on on 
the cooling have already been discussed in the literature 
long ago (Nomoto & Tsuruta 1987, Gudmundsson et al. 
1983). The models with GR effects included are marked as 
'GR' in Fig. || and drawn with thick lines. It is seen from 
those plots that the decay of the field is faster in NSs built 
on a soft EOS than in the case of a medium or stiff EOS 
but the decelerating GR effects are more pronounced. This 
is most remarkable for long living fields: the difference of 
the surface field strength for model 3, with the soft EOS, 
after 10"'^''yrs of evolution is about two orders of magni- 
tude. Comparing the field evolution for the soft and stiff 
EOS cases, while the final surface field strength for model 
3 in the stiff case is larger than those of the soft case by 
a factor of 300 when GR effects are neglected, that factor 
reduces to about 7 when relativistic effects are taken into 
account correctly. 
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In the late photon cooling era, when most of the ini- 
tial thermal energy of the NS has been irradiated away, 
the Joule heating by the decay of the crustal field com- 
pletely determines the cooling behavior. Notice, however, 
that most of the magnetic energy has been dissipated ear- 
lier when it had no effect on the cooling. The amount of 
heat released in the process of Joule heating is determined 
by the strength of the field at that time and its decay rate. 
Therefore, an initial field configuration as given by models 
of class 3 will result in a significant Joule heating while 
in models of class 1 the effect is practically nil. The field 
decay of the class 2 models yields some noticeable Joule 
heating only in case of the stiff EOS since its strength is 
large enough for that until about lO^yrs. For long peri- 
ods Joule heating is especially effective in NSs with a stiff 
EOS: their crust is quite thick, they keep a larger field and 
there is a lot of magnetic energy to dissipate. The effect of 
GR on the Joule heating is not spectacular but still signif- 
icant in the case of a soft EOS at the latest stages. For the 
different EOSs it can result in stronger or weaker heating: 
the rate of dissipation of magnetic energy is lower with 
GR but the magnetic energy to be dissipated is larger due 
to the previous slower evolution and the net effect can be 
an enhancement or a reduction of the heating. 

Notice, finally, that at these late ages when the thermal 
evolution is entirely controlled by the Joule heating the 
outer boundary condition, Eq. has no influence on the 
thermal evolution: this is due to the fact that the heating 
mechanism is temperature independent since a is totally 
controlled by impurity scattering. The star's luminosity 
is thus simply given by the total heating rate integrated 
over the whole star 



H 



qhe 



2* 



dV 



(31) 



since, given the low temperature at this times, qi, <C qii- 
This is very fortunate for our study since the Tf, — Tg rela- 
tionship' is poorly known at the low temperatures reached 
in this late phase of evolution. It would of course not be 
the case if the heating mechanism were temperature de- 
pendent. 

Notice finally that, in the case of magnetars, the Joule 
heating dominates the cooling from the very beginning 
and thus it has a strong effect onto the field evolution 
(Geppert, Page, Colpi & Zannias 1999). 

4-3. Magneto-thermal evolution within the 'fast' cooling 
scenario 

Enhanced neutrino emission is caused by the presence 
of exotic states of matter in the NS core. Kaon or pion 
condensation, hyperons, quarks as well as possible direct 
Urea processes enlarge the neutrino emissivity consider- 
ably. This results in a very fast cooling of the core, so 
that during the so called isothermalization phase the heat 
of the hotter crust is transported into the core. Depend- 



ing on the EOS and the assumptions about pairing that 
phase can last of the order of 1 to 100 yrs (see, e.g.. Page 
1998). During this phase the surface temperature drops 
very rapidly and after this the crustal temperature is so 
low that the conductivity is almost completely controlled 
by electron-impurity collisions. Thus, for a period of time 
given by the impurity decay time scale Tjmp, the field suf- 
fers practically no decay. For t > rimp the field decays 
according to a power law and, when it has diffused down 
to the crust-core boundary, the decay becomes exponen- 
tial. In class 3 models the power law like decay is missed 
because the field is already initially located close to the 
crust-core boundary so that it reaches this depth during 
Timp- Notice that Timp (x Q^^ and that it is also reduced 
with the softening of the EOS as an effect of the decreas- 
ing crustal thickness. Generally, the field decay is slower 
than in the 'standard' cooling scenario. 

Due to the accelerated cooling, the Joule heating, even 
for the model 1 and 2 fields, dominates the cooling ear- 
lier compared to the 'standard' cooling scenario. At ages 
above ~ lO"^ - 10^ yrs, however, the star's temperature is 
comparable, or even higher, to that predicted by the anal- 
ogous models within the 'standard' cooling scenario. In 
class 2 and particularly class 3 models with the Stiff EOS, 
the stars' temperatures in this range of ages are notice- 
ably higher, compared to the analogous 'standard' cooling 
cases, since the magnetic field is higher and thus the joule 
heating more efficient. After 10^°yrs, however, the differ- 
ence in the surface temperatures between the 'standard' 
and the accelerated cooling scenario vanishes. 

4.4- Rotational evolution 

Given the temporal evolution of the magnetic field for the 
various EOSs and cooling scenarios, the rotational period 
of the NSs can be estimated by integrating the equation 
which relates the loss of rotational energy to the radiation 
of electromagnetic energy by magneto-dipole radiation 



PP = 



27r^ R'^B^jt) 
~3 ^ ' 



(32) 



where B is surface magnetic field at the magnetic pole, and 
/ is the star's moment of inertia. A very accurately rela- 
tivistic expression for the moment of inertia, is given by 
/ = 0.21 £-2* MR^ (Ravenhall & Pethick 1994). However, 
this classical magneto-dipole radiation formula assumes 
flat space-time. Ideally we would like to have the exact 
GR version of this equation. To our knowledge, the GR 
version of it remains to be calculated and in the absence 
of the exact formula we shall only use an approximate ex- 
pression. We may recall that Eq. ^ implicitly considers 
the energy loss at the light cylinder and then extrapo- 
lates the field strength back to the star's surface assuming 
the flat space-time l/r^ radial dependence of the field. 
The full GR solution of a slowly rotating dipolar field, in 
vacuum, shows however (Muslimov & Tsygan 1990, 1992) 
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For each EOS, cooling scenario and initial field strength these three evolutionary tracks range from the most slowly 
decaying field (continuous curve) to the fastest decaying field (dashed curve) and thus encompass the whole range of 
predictions for the dynamical evolution of a neutron star with a crustal field. 



that the actual field at the star's magnetic pole is am- 
plified by a factor / = -3y^[y\n{l - y-^) + \{2 + y-^)], 
where y = R/Rs, compared to the flat space-time one. 
For instance, for our soft EOS star of 1.4 Mq , the re- 
sulting value of the amplification factor is / = 1.86 while 
for the stiff one, / = 1.27 (Page & Sarmiento 1996). In 



view of the Muslimov-Tsygan results, we will henceforth 
relate the surface magnetic field B at magnetic pole as 
we numerically calculated, it to the corresponding field B 
used in Eq. ^ for the rotational evolution hy B = f B. 
For each of our three EOS (and a fixed stellar mass of 1.4 
Mq ) and the two cooling scenarios, we specify the ini- 
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tial field strength and current depths po and finally the 
impurity content Qimp- The specification of those param- 
eters yields a unique evolution for the surface B-field and 
via Eq. ^ the corresponding P. The resulting magneto- 
rotational evolution of isolated NSs is presented in Fig. 0. 

We have also plotted, as dots, the PSRs data of the 
Princeton catalogue (Taylor, Manchester, & Lyne 1993): 
their B field is calculated with Eq. using P and P 
from that catalogue and the appropriate GR '/ factor' 
for each stellar model. The star's radius R and moment 
of inertia / are correspond to values of the correspond- 
ing models having mass of 1.4 Mq . The set of B-P dia- 
grams of Fig. 1^ is the arena where our theoretical models 
confront the reality. There is a major effect on the data 
interpretation arising from varying the EOS: we see that 
given the observed P-P data the "observed" surface B- 
field of pulsars differ almost one order of magnitude when 
we consider models where the equation of state is varied 
from the extreme soft one up to the extreme stiff one (see 
also UK97). This effect is a straightforward consequence 
of the fast growth of the magnetic field at small distance 
from the star which is moreover amplified by GR effects, 
combined with the much smaller stellar radius obtained by 
employing the soft EOS . It should emphasized however, 
that although a complete and accurate model of pulsar 
spin-down will probably alter the simple Eq. ^ (expected 
at least to modify the overall numerical coefficient) it cer- 
tainly will not change the radial dependence of the field 
in the near zone. 

A basic feature emerging from the Fig. (7), is the fol- 
lowing: a given model is viable if it manages to maintain 
a strong enough field for a long enough time such that 
its rotational period can increase up to values compatible 
with the bulk of the observed pulsar. It is also clear that if 
a 'standard' cooling scenario applies and irrespectively of 
the EOS, the field decays much faster than in the case of 
the accelerated cooling: this results in smaller saturation 
values of the rotational periods compared to the acceler- 
ated cooling model, with the same initial field structure 
and the same EOS. 

The evolutionary tracks offer the possibility to check 
the whether a given EOS can be naturally compatible with 
the observational constraints without extreme "fine tun- 
ing" of the initial parameters po and Bq, of the impu- 
rity parameter Qimp and its preference to standard ver- 
sus enhanced cooling. Thus for instance: assuming a soft 
EOS, a comparison of the evolutionary tracks with the 
observed pulsar population strongly favors the enhanced 
cooling scenario. Within the 'standard' cooling scenario, 
only the strongest initial fields ~ 10^^ G are acceptable. 
In this case, even Bq ^ 10^^ G requires very high po and 
very low Qimp, definitely making the 'standard' cooling 
scenario with a soft EOS very unappealing. Notice that 
the fast cooling scenario is very natural for a soft EOS 
since the central densities reached in NSs are above ten 
times nuclear matter density. Even if the pulsar fields (for 



a given P) were much lower than estimated we would still 
be in the same situation since the high field is required 
to be able to spin down the pulsars to the range of ob- 
served P's before the field decays too much. In the case of 
a medium EOS the requirements are much less stringent. 
Within the 'standard' cooling scenario, initial fields of the 
order of 10^^ G require high initial depths po and low pol- 
lution in order for the evolutionary tracks to reach the 
bulk of the pulsar population. Finally, in the case of a stiff 
EOS, the field decay is slow enough that the whole pulsar 
population is reachable with initial fields in the range 10^^ 
- 10^"^ G without basically any significant restrictions on po 
and Qimp in both the fast and the slow cooling scenarios. 

One cannot overemphasize that the above analysis is 
plagued by the intrinsic uncertainties of Eq. ^ and should 
be considered as only indicative. Moreover, it is worth 
stressing here that tracks starting with an initial field 
Bq — lO^^'G have to be considered with reserve since 
for such (and larger) field strengths the decay may not 
any longer described by a linear diffusion equation. The 
possible occurrence of a Hall-cascade (Goldreich & Reis- 
senegger 1992) would imply that the field evolution would 
deviate significantly from Eq. which is strictly valid in 
the limit that the magnetization parameter ujbt < 1 {ljb 
being the Larmor frequency and r the relaxation time of 
the electrons as given by Eq. ^0| ). 

5. Discussion and conclusions 

We studied in detail, and for a large variety of possible 
models, the magnetic, thermal and rotational evolution of 
isolated NSs, assuming that their magnetic fields, and the 
currents supporting them, are confined to the stellar crust. 
Our calculations take into account, for the first time, all 
mutual effects of the thermal and magnetic evolution self 
consistently in a wholly GR formalism. 

5.1. Comparison with the work of Urpin and Konenkov 

Urpin & Konenkov (1997, UK97 thereafter) have pre- 
sented the most detailed study of the evolution of crustal 
neutron star magnetic field while Miralles, Urpin & Ko- 
nenkov (1998) completed the former work by the inclu- 
sion of Joule heating in the thermal evolution. It should 
be mentioned, however, that neither of these works incor- 
porated GR effects on the i?-field evolution. 

Our field evolution results, without GR effects incorpo- 
rated, are very close to those of UK97 (their curves labeled 
3, dashed lines, in their figure 2 correspond to our models 
2). The main differences are seen in the early evolution of 
the field (our phase a) : this can be easily attributed to the 
fact that UK97 started the field evolution at an age of 100 
yrs while we started it immediately at the star's birth. 
Since UK97 used isothermal stars for their modeling of 
the B-field evolution they could not model the early hot, 
and non isothermal, phase properly. In our calculations. 
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this time difference allows the initial currents to rapidly 
relax from their initial distribution (Eq. which does not 
fulfill the outer condition) during a phase of low conduc- 
tivity while the models of UK97 must do it in conditions 
of higher conductivity, i.e., their relaxation is much slower. 
However, if one contemplates the scenario in which crustal 
magnetic fields are generated by a thcrmomagnetic insta- 
bility during this early hot phase (Blandford et al. 1983, 
Urpin, Levshakov & Yakovlev 1986, Wiebicke & Geppert 
1996) it is equally reasonable to start the B-field evolution 
by pure ohmic decay at the end of this phase as UK97 
did. The slight differences seen in the strength of the field 
at the onset of the impurity scattering dominated phase 
(phase b) can be easily explained as due to slight dif- 
ferences in the cooling histories, particularly differences 
in the time at which photon cooling takes over neutrino 
cooling (the knee at ages around 10^ - 10^ yrs) and also 
as a result of the differences in the previous phase a. 

Finally, the field evolution in the late phases, c and 
d, agree well with the results of UK97 once we take into 
account the small differences in the previous phases. 

5.2. General Relativistic effects 

While previous studies considered the influence of GR ef- 
fects into the star's structure and on its thermal evolutions 
(see e.g. Page 1989), their incorporation into the magnetic 
field evolution had not been properly accounted for. In the 
present paper we have been working with the GR version 
of the diffusion equation, Eq. |2^, accompanied by proper 
boundary conditions, Eq. derived in detail by GPZOO. 

The analysis of GPZOO clearly showed that most of the 
effect of GR on the field decay is due to the presence of 
the red-shift factor e* in Eq 



location of the currents and the impurity contents of the 
crust are significantly weakened when GR effect are taken 
into account. 

5.3. Comparison with the work of Sengupta 

As far as the previously mentioned work of Sengupta 
(1998) is concerned, although we are in qualitative agree- 
ment with his results, we find many quantitative as well as 
interpretational differences. For instance this author con- 
sidered a soft EOS (similar to the one we used) and field 
configurations with currents supposedly initially located 
at densities of 2 x 10^^ and 4 x 10^^ g cm~'^. However, we 
obtain much higher densities at the depths at which he 
locates these layers: using his values oi x = r/R ~ 0.979 
and 0.9834 we find densities of the order of 6 x 10^^ and 
2.6 X 10^^ g cm""^ respectively. Consequently, his field de- 
cay curves should be compared with our class 2 models: 
we see then a rough quantitative agreement with respect 
to the importance of the GR effects with the significant 
difference that the late exponential decay (phase d) is ab- 
sent in his GR models. As mentioned above, since the 
B-field evolutionary curves with and without GR effects 
should show approximately similar shapes we deduce that 
the different behavior of the B field at late times obtained 
by Sengupta, must be due to numerical inaccuracies com- 
bined with his neglect of appropriate boundary conditions 
at the stellar surface as well as his employment of the 
Schwarzschild geometry. 

5.4- Constraining the neutron star EOS and cooling 
history within the crustal magnetic field hypothesis 



20, This can be intuitively The analysis of the B — P diagrams of § 4.4 may be a 



understood by noticing that the red-shift factor relates the 
proper time in each layer inside the star, i.e., the physi- 
cal time in the layer where the currents are located and 
decaying, to the time of an observer who is observing the 
decay, at infinity (i.e., the coordinate time). As it can be 
easily seen from Eq. the effects of the red shift fac- 
tor on the field decay could be approximated by a flat 
spacetime diffusion equation running however on a slower 
time scale, given by some kind of averaged red-shift. This 
is clear from our Figures || and || where the GR curves 
have the same shape as the non-GR ones but shifted to 
the right. 

Comparing the field evolution when GR effects are in- 
cluded with the evolution when they are neglected, we 
see that we obtain quite larger fields already in phase a, 
and in the late decay (phase c), we obtain much larger 
fields, by up to two orders of magnitude for very compact 
stars. GR effects in a field evolution can be viewed as al- 
most equivalent to an evolution without GR effects but 
with currents initially located to higher densities and, in 
the late decay phase, with lower impurity content. This 
means that all constraints previously obtained about the 



tool to constrain the structure of matter at high density. 
In a similar analysis, UK97 concluded that a stiff EOS 
with 'standard' cooling is the most promising model for 
understanding the observed pulsar population properties 
and that a medium EOS requires currents located at rel- 
atively high densities and low impurity contents. These 
authors prefer the 'standard' coohng scenario on the basis 
that it implies an early decay of the field which may have 
an observational support coming the fact the young pul- 
sars with an associated supernova remnant have stronger 
magnetic field than the bulk of the population. Our re- 
sults, with GR included, show that a medium EOS with 
'standard' cooling is also compatible with the observed P 
and P with much weaker constraints on the initial dipole 
strength B^ and penetration density po and on the impu- 
rity content Qimp than when GR effects are neglected. A 
soft EOS with 'standard' cooling needs very special con- 
ditions to accommodate the observational data: large Bq, 
high po, and very low (3i,„p. 

But, notice that in case fast neutrino cooling is op- 
erating any EOS could be compatible with the observed 
pulsar population with only weak constraints on Bq, po, 
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and Qimp as we argued in § 4.4. However, this cooling 
scenario is more physical in the case of a soft EOS and 
probably incompatible with an EOS as stiff as our stiff 
EOS (for which the central density of a 1.4 Mq star does 
not even reach twice the nuclear density). In summary, it 
appears rather difficult within the crustal field hypothe- 
sis alone to draw any strong conclusion about the EOS of 
neutron stars and their cooling histories. 

5. 5. The EOS and cooling history of neutron stars and 
the crustal magnetic field hypothesis 

There exist independent arguments in favor of a not 
too soft neutron star EOS, none of them being, at the 
present time, compelling. Interpretation of kilohertz quasi- 
periodic oscillations (KHz QPO's) in several low mass X- 
ray binaries indicate these systems contain neutron stars 
with masses around 2 Mq (Kluzniak 1998) which a very 
soft EOS is not able to sustain. The strong gravitational 
light bending around a very compact star would make it 
almost impossible for such a star to show any modula- 
tion in its surface thermal emission, in contradistinction 
to some observations (Page & Sarmiento 1996). 

With regard to the thermal evolution of neutron stars, 
comparison of cooling models with current estimate of sur- 
face temperature of young neutron stars shows no clear ev- 
idence of occurrence of fast cooling (Ogelman 1995; Page 
1998). This means that if these neutron stars do contain 
some 'exotic' phase of matter, its strong neutrino emission 
must be quenched by pairing (Page & Applegate 1992; 
Page 1998) so that their thermal evolution is very close to 
the 'standard' one. 

These two lines of arguments are consistent with the 
conclusions arising from the crustal field hypothesis, in- 
cluding GR effects: the EOS of neutron stars is away from 
the very soft regime and their thermal evolution is close 
to the prediction of the 'standard' model. 

5.6. Joule heating and detectability of 
Old Isolated Neutron Stars 

Our study of the effect of the Joule heating produced by 
the decaying currents gives results very similar to the ones 
of Miralles et al (1998) and shows that the GR effects on 
the field evolution do not introduce any important change 
on the resulting late time thermal evolution. 

Another consequence of the GR effects on the B-field 
evolution is that we can predict significantly stronger field 
strength over a Hubble time. This implies that Old Iso- 
lated Neutron Stars will be spinning very slowly and are 
hence more likely to be able to accrete matter from the 
interstellar medium. As a consequence this significantly in- 
creases the chances of detecting them through their ther- 
mal radiation, either due to the Joule heating or to the 
accretion. 



6. Conclusions 

The effect of GR on both the thermal and magnetic evo- 
lution is to slow it down, mostly because the proper time 
inside the star runs more slowly than the observer's time. 
This effect is of course stronger the more compact the 
star is. On the other hand, with increasing compactness 
the field decay, as any diffusion process, is accelerated be- 
cause of the reduction of the length scale. The competition 
of these two opposite tendencies reduces the sensitivity of 
the field evolution to the softness or stiffness of the EOS. 
As a result it is difficult to draw conclusions about the 
nature of the dense matter EOS from magnetic field evo- 
lution studies alone. However when taking into account 
information from other approaches, as the cooling history, 
gravitational lensing, kHz QPO's, a consistent picture of 
neutron star structure and evolution is obtained, in which 
the EOS is not too soft and the cooling history is close to 
the so called 'standard' cooling scenario. 

Reversely, we can consider these results as an argu- 
ment in favor of the crustal magnetic field hypothesis. It 
is seen to be compatible with the observed distribution of 
pulsars in the P-P diagram without requiring fine tuning 
of the models' parameters. Moreover, these model parame- 
ters are consistent with values deduced from various other 
neutron star studies. 



Appendix A: Induction Equation and Joule 
Heating Rate 

In this appendix, we shall provide a few of the intermedi- 
ate calculations leading to the derivation of the relativistic 
expression for the Joule heating employed in the numeri- 
cal computations. We begin by recalling that the covariant 
form of Maxwell's equations are as follows (Misner et al. 
1973; Wald 1984): 



47r 



(A.l) 
(A.2) 



where Fap ~ —Fpa, Ja and V are the coordinate com- 
ponents of the Maxwell tensor, the conserved four cur- 
rent and the covariant derivative operator, respectively. 
We may first recall that once a solution Fa/3 of Eqs. [A.l| - 
|A.2| has been specified, an observer with four velocity 
jja ^ U"Ua = — 1 measures electric and magnetic fields 
{E , B) given respectively by: 



Fa — FapU^ , Ba — — -^^aP^^ F^sU^ 



(A.3) 



where eap^s stands for the four-dimensional Levi-Civita 
tensor density. It follows then easily from Eq. A.3 that an 
inversion yields: 



FafJ — UaFp - UpEa + eap-ysU^ 



(A.4) 
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Thorne et al. (1982, 1986) have introduced an elegant 
version of curved spacetime electrodynamics, the abso- 
lute space formulation, that is reminiscent of the famil- 
iar flat spacetime electrodynamics formulated in terms of 
the electric and magnetic fields. This can be accomplished 
by working directly with the physical frame components of 
the electric and magnetic fields {E , B), as determined by 
the static observers relative to their orthonormal frames. 
Recalling that such observers have a four velocity field 
= e*5^ (see Eq. 0) then it follows that Maxwell's 



Eqs. A.l- A.2 can be written in the following equivalent 
form: 



V-E = Anp , V-B = 

V>c{ZB) = ^ZJ+l^ 

c c at 



Vx{ZE) 



1 dB 

c~dt 



(A.5) 
(A.6) 

(A.7) 



where in this appendix Z stands for the lapse function re- 
lated to the red shift factor by Z = e* and the V operators 
are formed out of the following three metric: 



1 - 2Gm/c^r 



(A. 



by the rules of vector calculus as applied to the above 
orthogonal coordinate system. It follows then rather easily 
that in the absence of any convective motion Eqs. A.5-A.7, 
combined with Ohm's law J = aE, and within the MHD 
approximation (i.e., neglecting the displacement current in 
Eq. |A.6D , yield the following induction equation governing 
the evolution of the neutron star magnetic field: 



IdB 
c~dt 



Vx[^Vx(ZB)]=0, 



(A.9) 



The above form of the induction equation combined with 
Eq. A. 8 as well the line element Eq. |l| with the metric 
functions corresponding to a non singular perfect fluid so- 
lution of Einstein's equations yield Eq. of the main text 
(for a detailed derivation see GPZOO). 

We now consider the GR expression for the Joule heat- 
ing. To do so, we start from the covariant expression of the 
energy momentum tensor for an arbitrary electromagnetic 
field (Misner et al. 1973; Wald 1984) i.e.: 



1 

4tt 



(A.IO) 



Using the representation of Eq. A. 4, one can easily show 
that the electromagnetic energy density £ as seen by the 
static observers is given by: 



1 



— [E-E + B-B] =£ 
Btt 



(A.ll) 



where E and B stand for the physical components of 
the electric and magnetic fields, respectively, satisfying 



Eqs. A.5-A.7. Poynting's theorem now has the following 
form 



where 



2gS-jE 



S = —(E X B) 

An 

is the Poynting vector and 
\7Z 



(A.12) 



(A.13) 



(A.14) 



is the gravitational acceleration, dr — Z dt being the 
proper time interval. The second term in the r.h.s. of 
Eq. |Al2 is a purely relativistic effect which results from 
the inertia of (electromagnetic) energy (see, e .g., T horne 
et al 1986). The only dissipative term in Eq. A.12 is the 
last one, i.e., the Joule heating term, using Ohm's law: 

g(GR) ^ 21 ^^^^5^ 

By construction Q\^^^ represents energy per unit (proper) 
time and (proper) volume. 
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